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ABSTRACT 

We use iV-body simulations to study the infall of dark matter haloes onto rich clusters 
of galaxies. After identification of all cluster progenitors in the simulations, we select 
those haloes which accrete directly onto the main cluster progenitor. We construct the 
mass function of these merging satellites, and calculate the main orbital parameters 
for the accreted lumps. The average circularity of the orbits is e ~ 0.5, while either 
radial or almost circular orbits are equally avoided. More massive satellites move along 
slightly more eccentric orbits, with lower specific angular momentum and a smaller 
pericentre. We find that the infall of satellites onto the main cluster progenitor has a 
very anisotropic distribution. This anisotropy is to a large extent responsible for the 
shape and orientation of the final cluster and of its velocity ellipsoid. At the end of 
the simulations, the major axis of the cluster is aligned both with that of its velocity 
ellipsoid, and with the major axis of the ellipsoid defined by the satellite infall pattern, 
to as 30° on average. We also find that, in lower mass clusters, a higher fraction of 
the final virial mass is provided by small, dense satellites. These sink to the centre 
of the parent cluster and so enhance its central density. This mechanism is found 
to be partially responsible for the correlation between halo masses and characteristic 
overdensities, recently highlighted by Navarro, Frenk & White (1996). 
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£h ! 1 INTRODUCTION 



In cosmological models of gravitational instability, struc- 
tures form by the collapse of small perturbations of some ini- 
tial density field. The most successful flavour of these models 
is the hierarchical clustering picture, where matter clusters 
on small scales first, and structures form in a bottom-up 
fashion. In this scenario, clusters of galaxies represent the 
most recently assembled structures, and are therefore dy- 
namically young systems. For this reason, they are ideal 
candidates to study the link between the final structure and 
morphology of an object and its formation process. In this 
context, it is of great interest to investigate in some detail 
the accretion of matter within the proto-cluster. The study 
of the recursive merging of smaller size objects onto the main 
cluster progenitor is a powerful tool to unveil the cluster 
dynamical history, and can possibly help us to understand 
some details of its final structure. 

A detailed investigation of the cluster merging history 
usually requires either Monte-Carlo techniques or numer- 
ical simulations. In the former case, one uses an analyti- 
cal model for the clustering of structures, like the Press & 
Schechter (1974) (hereafter P&S ) formalism, and extensions 
(Bond et al. 1991; Bower 1991; Lacey & Cole 1993; Lacey 



& Cole 1994), to produce Monte Carlo realizations of the 
merging histories of dark matter haloes. These in turn are 
the starting point of a recipe for galaxy formation, whose 
predictions can be compared to the observed properties of 
real galaxies and clusters. Extensive work in this direction 
has been done in the past few years (Kauffmann & White 
1993; Cole et al. 1994; Heyl et al. 1995; Kauffmann 1995), in 
an attempt to interpret the various observations in a unified 
and global framework. While this work provides a plausible 
link between the properties of dark matter haloes and those 
of galaxies and clusters, it is necessary to complement the 
investigation with numerical simulations. These, although 
limited to a smaller number of applications, do not need 
as many simplifying assumptions on the details of structure 
formation. Therefore one can address issues not accessible 
to Monte Carlo techniques, as the dynamical equilibrium of 
clusters, the study of their density and velocity profiles, the 
characterization of the orbits of infalling haloes, the survival 
times of haloes after merging with the main cluster, and oth- 
ers. These results are also useful to improve the Monte Carlo 
methods, as they can tell which assumptions are more real- 
istic for the galaxy formation recipe. 

Although many numerical simulations of cluster forma- 
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tion have been performed in the past (White 1976; Quinn, 
Salmon & Zurek 1986; Navarro, Frenk & White 1995; Tor- 
men, Bouchet & White 1997) there is little work on the con- 
nection between the details of the cluster merging history 
and the final cluster configuration, in a cosmological context. 
One reason for this is that such analysis requires simulations 
with very high mass and force resolution, as only recently 
have become available. In the present paper we present the 
results of such a study. Using a sample of nine dark matter 
haloes of rich galaxy clusters, obtained from high resolu- 
tion TV-body simulations, we identify the merging history of 
all cluster progenitors, and extract from that a population 
of satellites, i.e. progenitor haloes accreting directly onto 
the most massive cluster progenitor at any time. We study 
this population and characterize its dynamical properties. 
We then use the information on the infalling satellites as a 
key to interpret the cluster formation history, and its final 
shape and structure. In a forthcoming paper (Tormen 1997, 
in preparation) we shall make a detailed comparison of the 
merging history of the clusters with the predictions of the 
extended P&S formalism. 

The outline of the paper is as follows. Section 2 briefly 
presents the simulations, and describes the method used to 
define their halo population. In Section 3 we study the clus- 
ter satellites. We look at the way they merge with the main 
halo, and measure their mass distribution. Then, by moving 
the satellites in the spherically averaged, static potential of 
the main halo, as measured at each time output, we inte- 
grate their orbits and derive the main orbital parameters. 
In Section 4 we study the connection between shapes and 
orientations of clusters and the infall pattern of their satel- 
lites, and investigate the dynamical origin of universal dark 
matter density profiles (Navarro, Frenk & White 1996). Fi- 
nally, in Section 5 we summarize the paper and present the 
main conclusions. 



2.2 Identification of dark matter lumps 

The relatively low number of particles of each simulation 
(N ~ 60000 on average) allowed us to use the potential 
energy of particles to identify the location of dark matter 
lumps (Efstathiou et al. 1988). We defined lumps by an over- 
density criterion, and included all particles within a sphere 
of mean overdensity 5 V = 178, centred on the particle with 
lowest potential energy. The value 8 V corresponds to the 
virial overdensity in the model of a spherical top-hat col- 
lapse. We shall call virial radius r v the corresponding radius 
of the lump. We identified all lumps with at least five parti- 
cles within r v . 

In order to determine the minimum size of lumps to 
use in our analysis, we performed the following test. We se- 
lected, for each output of each simulation, all lumps formed 
by n particles, and looked at the fate of those particles in the 
next time output. The particles of a lump will now be part 
of other lumps. We took note of the largest fraction of lump 
relic that was found in a single halo, and looked at the dis- 
tribution of such quantity. If all, or most, of a lump is found 
in another halo on the next time output, the lump is a dy- 
namically consistent object. We performed the analysis for 
several values of n, from n = 5 up, and found convergence in 
the resulting main relic distribution for lumps with at least 
8 particles. For such lumps, in 95 per cent of the cases the 
largest relic carried at least half of the original mass, and in 
70 per cent of the cases it carried more than 80 per cent. On 
the other hand, almost 20 per cent of lumps with 5 particles 
were completely dissolved in the field by the next output. 
Therefore we chose to limit our analysis to lumps with at 
least eight particles within their virial radius. All remaining 
particles, including those in smaller lumps, will be generally 
referred to as field particles. Since we limit our analysis to 
lumps which are progenitors of the nine simulated clusters, 
these lumps are only made up with high resolution particles. 



2 METHOD 

2.1 Simulation setup 

For the present analysis we use the iV-body simulations pre- 
sented in Tormen et al. (1997), where complete details on 
the simulations may be found. In summary, our sample con- 
sists of nine dark matter haloes of galaxy clusters, with an 
average mass of 1.13 x 10 15 M Q , resolved by « 20000 dark 
matter particles each within the virial radius, with an effec- 
tive force resolution of ~ 25 kpc. The haloes come from an 
Einstein-de Sitter universe with scale-free power spectrum 
of fluctuations P(k) oc k n , and a spectral index n = — 1, 
the appropriate value to mimic a standard cold dark matter 
spectrum on scales relevant to cluster formation. The sim- 
ulations have a Hubble parameter Ho =50 km s _1 Mpc -1 
and are normalized to match the observed local abundance 
of galaxy clusters (White et al. 1993). Each simulation was 
performed using two sets of particles, a first high resolution 
set to actually form the cluster, and a low resolution set of 
massive background particles to model the large-scale tidal 
field of the simulation out to a scale L — 150 Mpc. For each 
of the nine simulations we shall use ten outputs, ranging 
from 2 ~ 6 to z — 0, for a total of 90 outputs. 



3 MERGING OF SATELLITES 

Let us call, at any given redshift, progenitors of a cluster all 
haloes containing at least one particle that by z — will 
be part of the cluster, and also all field particles which will 
end up in the final cluster. At any redshift, we define the 
main (or largest) progenitor of a cluster as the progenitor 
halo containing the largest number of particles which will 
be part of the final cluster. In all cases this coincided with 
the most massive progenitor halo, although, with the present 
definition, the main progenitor in a simulation is not neces- 
sarily the most massive one. To help clarifying the meaning 
of these terms, one can refer to the schematic merging his- 
tory of halos given by the merging tree in Fig. 6 of Lacey & 
Cole (1993). 

We termed satellites those progenitor lumps which were 
directly accreted onto the main cluster progenitor. Opera- 
tionally, we started from the list of progenitors at a given 
output of the simulation, and selected those haloes for which 
at least one particle was found in the largest cluster progen- 
itor in the next time output of the simulation. We shall 
characterize these satellites by studying their mass distribu- 
tion, and by calculating their orbital parameters. We shall 
call identification time Ud of a satellite the last output time 
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Figure 1. Number weighted distribution of the fraction of satel- 
lite mass which is captured by the main cluster progenitor. Each 
curve refers to a different range of particle number for the satel- 
lites: m v < 20 particles (solid), 20 < m v < 100 particles (dotted) 
and m v > 100 particles (dashed). The two extreme situations 
are most common: satellites stripped by a small amount of their 
mass, with m captured /m v < 0.2, and satellites captured (almost) 



entirely, with m captv 
rare. 



ed/ m v > 0.8. Intermediate cases are very 



before the satellite accretes onto the main cluster progeni- 
tor. Since we consider ten outputs for each simulation, there 
are nine outputs (all but the last) to identify satellites. 

In the hierarchical clustering picture of structure forma- 
tion, matter clusters on small scales first, and dark matter 
haloes of a given size are formed by the assembly of smaller 
haloes which were formed earlier. In this idealized model, 
when a satellite lump is accreted onto a more massive halo, 
all of the satellite mass is incorporated in the bigger object. 
Although this description is fairly close to what we observe 
in iV-body simulations of hierarchical clustering, the actual 
collapse of structures is more complicated than its idealized 
version; in particular, satellites may be tidally stripped by 
the accreting halo, so that only a fraction of the satellite 
mass becomes part of the accreting system. This is clearly 
illustrated in Fig. [l], where we show the number distribution 
of the fraction f c of satellite mass captured by the main clus- 
ter progenitor, for all time outputs and all clusters together. 
Each curve refers to satellites with different number of par- 
ticles. We observe that most of the satellites transfer either 
very little or most of their mass to the main cluster pro- 
genitor, whereas less than 20 per cent of the satellites have 
fc 6 [0.2, 0.8]. This trend is more evident for more massive 
satellites. While this result may be partially due to a few 
satellites which, by chance, happen to cross the virial radius 
of the main progenitor at the time of a simulation output, 
it is more likely that the distribution reflects different phys- 
ical situations. We shall address the issue at the end of the 
section on the orbital parameters. 
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Figure 2. Mass function of satellites. The histogram shows the 
fraction of cluster mass accreted in satellites of mass m v /M v (to), 
averaged over all clusters in our sample. Only haloes accreted af- 
ter the cluster formation redshift zt (the redshift when the most 
massive progenitor reaches half of the final cluster mass), and 
captured by at least 50 per cent of their virial mass, are consid- 
ered. The mass of satellites is given in units of the final cluster 
mass; error bars are Poissonian. The solid curve is a prediction 
from Monte Carlo merging trees, taken from Lacey & Cole 1993 
(their Fig. 12). 



3.1 Mass function of satellites 

It is interesting to look at the spectrum of masses that 
clumped together and formed the clusters in our sample. 
This will tell us, for example, how much of a cluster's mass 
was made up by violent encounters with massive satellites, 
as opposed to smooth accretion of small sized systems. We 
want to select from the satellite list those who significantly 
contribute with their mass to the formation of the cluster. 
We do this by eliminating all satellites captured by less than 
50 per cent in mass, that is with f c < 0.5. This choice seems 
reasonable, as Fig. |l] shows that there is a clear separation 
between satellites barely stripped and satellites really cap- 
tured. Moreover, with this distinction, a satellite excluded 
from the selection in one output may be included in the 
next, if it is really captured by the main halo. We shall call 
merging time of the satellite onto the cluster the time when 
the satellite first crosses the virial radius of the main cluster 
progenitor. Since the separation between consecutive out- 
puts of our simulations is At ~ 1.6 Gyr, we can only bracket 
the estimated merging time of a satellite between the iden- 
tification time tid and the next output of the simulation at 
U d + At. 

Fig. jsj shows the mass distribution of satellites captured 
by at least 50 per cent in mass; the abscissa is the virial (i.e. 
total) mass of the satellite normalized to the final mass of 
the cluster. The smooth curve is taken from Lacey & Cole 
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Figure 3. The same distribution shown in Fig. ^] is shown here 
with satellite masses normalized to the mass of the main cluster 
progenitor at the time of capture. Both the cumulative (solid line) 
and the differential (thin dotted line, with Poissonian error bars) 
are drawn. These curves show the relative occurrence of merging 
events with satellites of a given mass ratio m v (z)/M v (z). 



1993, and is the expected mass distribution of satellites ob- 
tained from Monte-Carlo merging histories of haloes. The 
prediction does not depend on the final cluster mass if one 
considers only the accretion events happening after the for- 
mation redshift Zf, defined as the redshift when the most 
massive cluster progenitor first reaches at least half of the 
final cluster mass (Lacey & Cole 1993, 1994). In our case, 
we only considered satellites with Zid < Z/. To give an idea 
of our statistics, there are roughly 1400 satellites merging 
with the clusters, of which almost 600 identified after Zf. 
The histogram is normalized to the fraction of cluster mass 
provided by the satellites with f c > 0.5. This is done to keep 
account of the matter accreted onto the main cluster pro- 
genitor in both field particles and satellites with f c < 0.5. 
The agreement between model and simulation data is very 
good, as the curve matches the histogram over nearly three 
orders of magnitude in mass. 

In Fig. 0, the a;— axis has masses normalized to the final 
mass of the clusters, for the purpose of comparing it with 
the Monte-Carlo prediction. However, it is more instructive 
to look at the same distribution normalizing masses to that 
of the main cluster progenitor at the identification time of 
the satellite, tid- This is done in Fig. 0, where the differential 
(dotted histogram) and cumulative (solid histogram) mass 
functions of the accreted satellites are plotted. Now we use 
all satellites identified at any tid, but always with f c > 0.5. 
The (Poissonian) error bars refer to the differential distribu- 
tion. 

In the mass range log x <; —2 the distribution follows 
a power law with slope 0.5, the same slope of the distribu- 



tion in Fig. [2|. Normalization is the same described above, 
showing that, on average, one third of the cluster mass was 
accreted either in field particles or in satellites captured by 
less than half of their mass. The next largest contribution, 
about 20 per cent of the final cluster mass, comes from en- 
counters with satellites with mass ratio m v /M v > 0.5. On 

the whole, roughly 40 45 per cent of the total cluster mass 

comes from encounters with satellites with m v /M v > 0.1. 
These figures are quite robust, as they remain essentially 
the same if we use in the histograms all satellites (regardless 
to the captured fraction of their mass) or if we use for the 
satellite mass the actual captured mass instead of the virial 
mass. How frequent are encounters with massive satellites? 
Let us consider only satellites with f c > 0.5. The average 
number per cluster of encounters with satellites with mass 
ratio > 0.1 (> 0.2; > 0.5 respectively) is 2.7 ±0.5 (1.5 ±0.4, 
0.1 ± 0.1) from the time tj when a cluster reached half of 
its final mass (zf m 0.5 on average) to the present time, and 
13.7±1.4 (7.4±1.0, 2.3±0.5) over a Hubble time. Error bars 
are la estimates from bootstrap resampling over the clus- 
ters. It has been shown (Tormen et al. 1997) that encounters 
with satellites with mass ratio larger than roughly 0.2 can 
perturb the rms velocity of the cluster by as much as 30 per 
cent within the virial radius. According to the present re- 
sult, a rich cluster should undergo such encounters between 
one and two times since zj. 

3.2 Orbital parameters 

In this section we study the orbits of the satellites which 
accrete directly on the main cluster progenitor. Since we 
want to investigate how satellites fall onto the cluster, we 
are interested in the initial orbital parameters, which define 
the orbits prior to the actual merging. Of course, after a 
satellite crosses the virial radius of the main cluster, differ- 
ent dynamical processes (tidal stripping, dynamical friction, 
encounters with other halos, etc..) become more relevant, 
and its orbit will be modified. However, the initial orbits 
are still very interesting as they trace the accretion process. 
Moreover, they can be used as an input for a more realis- 
tic modelling of the dynamical processes listed above, for 
example in semianalitic models of galaxy formation. 

In order to integrate the orbits of satellites we used 
the static spherical potential approximation (e.g. Binney & 
Tremaine 1987, their Eq.(3.13)), and moved the satellites 
in the potential of the main cluster progenitor as measured 
at time tid- The radial potential profile <E>(r) was obtained 
by averaging the actual potential of all the particles in each 
output in spherical shells of equal logarithmic width 0.05, 
centred on the main cluster progenitor. A further Gaussian 
filter, of logarithmic width 0.2, was applied to the profiles 
to make them smoother and so avoid multiple solutions in 
the orbital parameters. Although we calculated the orbital 
parameters for all satellites, here we shall only show results 
for satellites with f c > 0.5. 

The approximations made will lead sometimes to unre- 
alistic solutions, especially (a) for massive satellites, whose 
contribution to the total potential is not negligible, and 
(b) at early times, when the main cluster progenitor grows 
faster, and so the static approximation is poorer, and when 
there is more chance of having big neighbouring haloes that 
invalidate the spherical approximation. In order to mini- 
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mize these problems, in this Section we limit our analysis 
to outputs where the main cluster progenitor has at least 
half of the cluster final mass, that is to redshifts za < Zf. 
With this choice, the mass of satellites, relative to that of 
the main cluster, is kept small (in our case m v /M v < 0.25 
for more than 99 per cent of the satellites), and the static 
model is closer to reality, since the potential well of the clus- 
ter has already been built up. This selection will somehow 
bias our statistics, because we automatically exclude mas- 
sive (m„/M„ <; 0.25) satellites from the analysis. However, 
we prefer this choice to the alternative of including unreli- 
able satellites in the results. 

In order to compare results from different outputs and 
from different clusters, we need to rescale them to uniform 
units. In particular, lengths should be expressed in units of 
the virial radius of the main cluster progenitor. Now, even in 
a static potential, the virial radius R v of the main halo will 
grow with time, since the mean background density of the 
universe is decreasing like (1 + z) 3 . As a consequence, from 
the time tj<j when we identify a satellite, and the later times 
t v , t p when the satellite crosses the virial radius and reaches 
the orbit pericentre, R v will have grown, sometimes by a not 
negligible amount. This growth must be taken into account 
since our satellites are identified at different distances from 
the main progenitor, and so cross its virial radius at different 
times. For every satellite, we calculate t v and R v (t v ) moving 
the satellite in the static potential of the cluster, and moving 
the cluster virial radius outward until the satellite crosses it. 

A proof that the assumed model is a better approxima- 
tion to the real clusters at late than at early times comes 
from the distribution of orbit pericentres: the most obvious 
unphysical result is a satellite not merging at all with the 
main cluster, that is, with a pericentre r p > R v . This actu- 
ally happens for 14 per cent of the satellites with zm > Zf, 
but only for 2 per cent of the satellites with Zid < Zf. 

We must also be careful that the satellites are not too 
far from the main cluster progenitor at their identification 
time. In fact, at large distances from the cluster, orbits are 
more easily influenced by neighbouring massive objects, and 
the spherical average of the potential is no more a good de- 
scription of the true field. Fortunately, in practice this is not 
a problem, since all selected satellites have nd < 3R v (tid), 
and 90 per cent have nd < 2R v (Ud). 

What is the shape of the orbits? How close do the satel- 
lites come to the centre of the main cluster? How much an- 
gular momentum do they carry? These are some questions 
we want to answer in this Section. Fig. [l| shows a scatter 
plot of the orbital circularity e versus the pericentre of the 
satellite's orbit (the latter rescaled to the virial radius of 
the main cluster progenitor, R v (t v ), at the time of crossing 
it). The circularity e = J/J C {E) (e.g. Lacey & Cole 1993) is 
the ratio of the angular momentum J of the orbit to that, 
J C (E), of a circular orbit with the same energy E. Radial 
orbits have e — 0, while e = 1 for circular ones. The figure 
shows that more eccentric orbits come closer to the cluster's 
centre, as naturally expected. The mean pericentric distance 
is r p /R v (t v ) = 0.38±0.26 (1 a of the distribution). The dis- 
tribution for e has an average value of 0.53 ± 0.23. Almost 
radial or almost circular orbits are much less likely than in- 
termediate ones. The mean amplitude of the 3-dimensional 
velocity of satellites as they cross R v is v aa t(Rv) = 1405±234 
km s _1 , with 5, 25, 50, 75 and 95 percentiles at 990, 1260, 
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Figure 4. Scatter plot of the orbital pericentre r p /R v and the 
orbital circularity e. The two histograms show the corresponding 
projected distributions (solid line for r p /R v and dashed line for 

1510, 1760 and 2200 km s _1 respectively. In terms of the cir- 
cular velocity of the main halo at R v (t v ), the same velocity 
has 5, 25, 50, 75 and 95 percentiles at 0.77, 0.99, 1.12, 1.21 
and 1.38 respectively. 

We are now in the position to check if satellites only 
partially captured and satellites completely merged have 
different dynamical histories. One possibility is that the 
formers are on more circular orbits than the latters, so 
that only the satellite outskirts are tidally stripped by 
the main progenitor, and the rest of the satellite escapes. 
At the other extreme, satellites in more eccentric or ra- 
dial orbits would be completely swallowed. For satellites in 
these two extreme situations we found average circularities 
e(f c < 0.2) = 0.55 ± 0.01 and e(/ c > 0.8) = 0.527 ± 0.007 
(la of the mean, now considering satellites at all ta). There- 
fore the difference is not significant enough to explain the 
effect of Fig. |lj Another possible explanation is that satel- 
lites with f c <C 1 have already had a first encounter with 
the cluster, but they had enough kinetic energy to emerge 
from it once, leaving behind only part of their mass. The av- 
erage 3-dimensional velocities when crossing R v are, for the 
same subsets, v sa t{Rv\ fc < 0.2) = 1370 ± 23 km s" 1 and 
v 3a t{Rv\.fc > 0.8) = 1239±12kms" 1 (la of the mean). Here 
the difference is more significant: indeed satellites captured 
by less than 20 per cent in mass have a higher speed, at a 3cr 
confidence level, than those completely captured. Perhaps a 
combination of this effect and of the former is responsible 
for the shape of the distribution shown in Fig. |l[ 

To see if there is any dependence of these and the other 
orbital parameters on the mass of the satellite (relative to 
the main cluster progenitor), we binned the data accord- 
ing to the relative mass of the satellite, m v /M v , where M v 
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Figure 5. Dependence of some orbital parameter on the satellite mass. Each panel shows the cumulative distibution of a parameter, 
in three mass ranges: m v /M v > 0.05 (solid lines); 0.005 < m v /M v < 0.05 (dotted lines) and m v /M v < 0.005 (dashed lines). Top 
Left: orbital pericentre r p /R v ; Top Right: amplitude of the angular momentum, J/J v \ Bottom Left: orbit circularity e; Bottom Right: 
radial period T r /Gyr. In each panel, the numbers shown are the result of a KS test between every pair of distribution: pi2 refers to the 
solid-dotted pair; P13 to the solid-dashed, and P23 to the dotted-dashed. 



is calculated at the time t v when the satellites cross R v . 
The results are shown in Fig. ^, where we plot the cumu- 
lative distributions of r p /R v , of e, of the amplitude of the 
angular momentum of the orbit relative to that of a circu- 
lar orbit at the virial radius, J/J v , and of the radial pe- 
riod T r of the orbit. Data are binned in mass, in the three 
intervals m v /M v < 0.005, 0.005 < m v /M v < 0.05, and 
m v /M v > 0.05. The figure shows some trend for more mas- 
sive satellites to be on slightly more eccentric orbits, which 
penetrate deeper in the cluster, with lower (specific) angular 
momentum J/ J v and with a shorter orbital radial period T r . 
We estimated the significance of these trends by a KS test 
performed on each pair of distributions. The probability of 
drawing such distributions by chance, from an identical un- 
derlying population, is given in each panel of the figure. This 
test shows that the trend with mass is statistically more sig- 
nificant for the radial period T r and the amplitude of angular 
momentum J/ J v , and less for the pericentric distance r p /R v 
and the shape e of the orbits. In the next Section we briefly 
discuss a possible explanation for these results, although we 
stress that these trends with mass are rather weak. An even 
weaker trend was also found for the kinetic energy of satel- 
lites as they cross R v : more massive satellites approach the 
main halo at slightly lower speed than less massive ones. 
The effect is small, less than 10 per cent for the velocities, 
and is possibly due to dynamical friction. We finally looked 



for dependence of these parameters on the mass of the final 
cluster, or on the time of identification tid, but did not find 
any significant trend in the results. 



4 DISCUSSION 

4.1 Collapse anisotropy and shape analysis 

Cosmological iV-body simulations have shown that, in most 
scenarios, the gravitational collapse of matter happens first 
along pancakes, or sheet-like structures. Then filaments 
form, possibly at the intersection of different pancakes. Mat- 
ter in these filaments collapses into dark matter haloes, and 
these flow along the filaments towards the potential min- 
ima, where they form galaxy clusters, perhaps at filaments' 
intersections. The collapse is therefore a highly anisotropic 
process. Recent observations of the distribution of substruc- 
ture in X-ray clusters (West, Jones & Forman 1995) seem 
to indicate that anisotropic collapse is common also in real 
clusters. West et al. (1995) suggest that this anisotropic for- 
mation may be responsible for the shape and orientation of 
clusters with the surrounding large-scale structure, on scales 
~ 10/i -1 Mpc, which is observed in different samples (e.g. 
Binngeli 1982; West 1989; Plionis 1994). A similar alignment 
is also found in cosmological iV-body simulations (West, Vil- 
lumsen & Dekel, 1991; Splinter et al, 1997). The link between 
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cluster shapes and merging histories has been addressed in 
the past (van Haarlem & van de Weygaert 1993). It is also 
recognized that some alignment exists between the cluster 
body and its velocity ellipsoid (e.g. Warren et al. 1992). 

The high resolution of our simulations, together with 
the knowledge of their complete merging history, allows us 
to address the issue in a more direct and quantitative way 
than in previous work. We proceed as follows. For each clus- 
ter of our sample, we first found the principal axes of the 
cluster mass distribution, and of its velocity ellipsoid, by di- 
agonalizing the mass tensor and the velocity dispersion 
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where k labels particles and i, j are the Cartesian compo- 
nents of the vectors. Since all particles are identical, their 
mass does not appear in the expression of Uj. The vector x 
identifies the cluster centre, and the mean velocity v is cal- 
culated using all the particles inside 0.5R V to limit the bias 
from infalling substructure. The sums are extended to all 
particles within the virial radius. Although this choice has 
been shown to systematically underestimate the true axial 
ratios (Warren et al. 1992), this is not relevant for the main 
purpose of the present analysis. 

To characterize the infall pattern we similarly defined 
the mass tensor of the satellite distribution: 

lij = £k,iXk,jiri k>v . (3) 

k 

Here x is the unit vector pointing in the direction (0, (j>) of 
approach of the satellite towards the main cluster progeni- 
tor. This direction is measured at time tid, just before the 
merging of the two. Each contribution to iy is now weighed 
by the satellite's virial mass m v . The reference system is 
centred on the main cluster progenitor. 

We then compared the shapes of the three ellipsoids 
and the relative orientations of the three major principal 
axes. Shapes can be classified using the triaxiality parameter 
(Franx, Illingworth & de Zeeuw 1991): 



(a 2 - b 2 ) 



(4) 



(a 2 - c 2 ) ' 

where a > b > c are the three principal axes of the el- 
lipsoid. We shall call oblate ellipsoids with < T < 1/3, 
triaxial those with 1/3 < T < 2/3, and prolate those with 
2/3 < T < 1 (Warren et al. 1992), and will label T Cl T v 
and Ts the triaxiality parameters referring respectively to 
the cluster shape, to its velocity dispersion ellipsoid and 
to its satellite distribution as defined above. We quantify 
the alignment in the orientation of two ellipsoids E\ and 
E2 using the cosine of the (3 dimensional) misalignment 
angle between the two major axes a\ and 02, defined by 
cos (A) = cos 61 cos 82 +sin 81 sin 82 cos(</>2 — (j>i). This is uni- 
formly distributed in [0, 1] for a random distribution of an- 
gles. Values larger than 0.5 (corresponding to misalignment 
angles smaller than 60°) indicate a correlation in orienta- 
tion. 

In order to study the time evolution of shapes and orien- 
tations, we calculated them at two different times: at redshift 
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Figure 7. Orientation correlations between the major axes of the 
cluster shape (C), its velocity ellipsoid (V) and the distribution of 
its satellites (S). Each panel shows the distribution of the cosine 
of the misalignment angle for a pair of axes, and is labelled with 
the mean and la of the mean. In the relevant cases we also give 
the probability p of obtaining the measured mean alignment from 
randomly oriented vectors. 



z = and at redshift z$ when the main cluster progenitor 
first assembled half of the final cluster mass. We only con- 
sider satellites identified at z < zj. We do so because for 
2 > Zf the most massive progenitor may be on different 
branches of the merging tree at different redshifts. In such a 
case one cannot relate the satellite infall to the halo shape. 
Moreover, at early times satellites are identified at larger dis- 
tances from the main halo, in units of its virial radius, and 
the direction of infall is less well determined due to possible 
perturbations on the orbits. 

The results of our analysis are presented in Figs [i] and Fl 
and in Table ^. They refer to satellites with f c > 0.2. Fig. g 
correlates, for each cluster of our sample, the triaxiality pa- 
rameter of the cluster shape, that of its velocity dispersion 
and those of the satellite distribution. Fig. ^ does the same 
for the cosine of the misalignment angle cos(A). Table [l] 
gives the ratio of major to minor axis for each considered 
ellipsoid. 

From the figures and the table we can make the follow- 
ing comments. The cluster shape and velocity ellipsoids are 
preferentially prolate or triaxial, and somewhat correlated 
in shape, at early times (Fig. ^a). As evolution proceeds, 
shapes do not significantly evolve. Although the correlation 
improves from Fig. |^a to Fig. ^b, this change is not statis- 
tically significant according to a standard KS test. Table [l] 
shows that the elongation of cluster shapes and of their ve- 
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Figure 6. Shape Correlations between the cluster (C), its velocity ellipsoid (V) and the distribution of its satellites (S). Each number 
labels a cluster. Panel (a) compares the shape of the cluster and that of its velocity ellipsoid at redshift Zf when the main cluster 
progenitor first assembled half of the final cluster mass. For our simulations Zf fa 0.5. Panel (b) does the same at z = 0. Panels (c) and 
(d) correlates the spatial and velocity structure of the cluster at z = with the distribution of infalling satellites for z < Zf . The grid in 
each panel separates region corresponding to different shapes: oblate (T < 1/3), triaxial (1/3 < T < 2/3) and prolate (T > 2/3). The 
error bars indicate the mean ±lcr of the mean for each distribution. 



Table 1. Ratio a/c of longest to shortest axis for the cluster shape 
ellipsoid (C), for its velocity ellipsoid (V) and for the ellipsoid 
defined by the infall pattern of satellites (S). For each ratio, the 
mean ±1<t of the mean is given. 

Label a/c a/c 

tf t 

C 1.57 ±0.08 1.60 ±0.06 
V 1.45 ±0.06 1.32 ±0.05 
S - 3.16 ±0.30 



locity ellipsoids do not appreciably change in time. A strong 
shape-velocity correlation is shown in the axes orientation, 
both at early and late times (Fig. 0a,b). By the present 
time the mean misalignment between the cluster shape and 
its velocity ellipsoid is of the order of 30°, with a very high 
statistical significance. 

On the other hand, both the cluster and its velocity 
ellipsoid change orientation in time, so that there is no very 



significant correlation between the cluster orientation at tf 
and that at to (Fig. file), and between the corresponding 
velocity ellipsoids (Fig. |7|i). 

The final cluster shape is strongly aligned (« 30° on av- 
erage) with the infall pattern of satellites accreted after tf 
(Fig. ^e). The shapes also are similar (Fig. |^c). A similar cor- 
relation is observed for the relative alignment of the velocity 
ellipsoid and of the satellite pattern (Fig. ^|f). The shapes 
are correlated in that they tend to avoid oblate configura- 
tions (Fig. The satellite pattern is prolate or triaxial 
in eight clusters out of nine. The axial ratio a/c in Table [I] 
indicates a very anisotropic infall. 

The same analysis on the orientation of the minor and 
intermediate axes of the ellipsoids gives a picture consistent 
with these results. In particular, it appears that the major 
and intermediate axes of an ellipsoid may sometimes switch 
during cluster evolution, causing anticorrelations similar to 
those observed in Fig. ^c,d. The minor axis instead remains 
rather well defined. 

To summarize, the present results show that the 
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anisotropic infall of matter on the forming cluster is to a 
large extent responsible both for the final cluster shape and 
for its velocity structure. Both are strongly correlated, in 
shape and orientation, with each other and with the infall 
pattern of merging satellites. In their analysis, van Haarlem 
& van de Weygaert (1993) found that the cluster orienta- 
tion is mainly determined by the last infalling subcluster. 
Although our results are in general agreement with theirs in 
highlighting the origin of such orientation, the present anal- 
ysis suggests that cluster alignments are stable enough to be 
observed for a long time, as the correlations we found reflect 
the whole cluster history after its formation at Zf ~ 0.5. 

Even if we have not directly checked the alignment 
of our clusters with the neighbouring large-scale structure, 
our analysis strongly suggest this picture, since the infalling 
satellites are tracers of the matter distribution around the 
clusters. Therefore our results are fully consistent with the 
alignments observed between clusters and the surrounding 
structure on scales ~ 10/i -1 Mpc, and favour an explanation 
of these in terms of structure accretion along preferred di- 
rections, e.g. along filaments. The shape and orientation of 
clusters would then be mainly determined by the large-scale 
initial density field. 

This conclusion may also explain the trend in the orbital 
parameters shown in Fig. ||. Massive satellites are more likely 
to reside in dense filaments than less massive ones. There- 
fore, the formers should approach the main cluster progeni- 
tor along more eccentric orbits than the latters, with lower 
specific angular momentum and a closer pericentric distance, 
as indeed found. 

4.2 Universal density profiles 

Navarro et al. (1996) have recently found that dark haloes 
from their cold dark matter simulations exhibit remarkable 
similarity in shape, over a wide range of masses. After rescal- 
ing lengths in units of the halo virial radius R v , they could 
fit the dark matter density profile of objects from the size 
of dwarf spheroidals to that of galaxy clusters using the one 
parameter fit 

P( r ) = 8- 
p b (r/r s )(l + r/r s ) 2 ' 1 ' 

The free parameter is the scale radius r s , which in turn de- 
termines the characteristic halo overdensity S s , defined by 
Equation (JH|). Navarro et al. (1996) found that, expressing 
halo masses in terms of the characteristic non-linear mass 
M* corresponding to a linear overdensity of order unity, 
more massive haloes have larger scale radii, i.e. larger values 
of r s f Rv, or equivalently have lower characteristic overden- 
sity S s - Further work (Cole & Lacey 1996; Tormen et al. 
1997, Navarro, Frenk & White 1996b) have confirmed this 
trend and extended the result to different cosmological mod- 
els. This behaviour is in fact consistent with the expectations 
of hierarchical clustering, where structures form in a bottom- 
up fashion. Since at earlier times the universe was denser, 
one expects that less massive haloes, which formed earlier, 
have higher characteristic density S s , while more massive 
haloes, which form later from the merging of less massive 
ones, have lower values of S s . Navarro et al. (1997) link the 
characteristic density of a halo to the epoch when its pro- 
genitors first collapsed. They find a good agreement between 



their numerical results and the analytical predictions of the 
extended P&S formalism. Note however that Navarro et al. 
define the formation redshift of a halo in a way different 
from ours. 

A physical mechanism for the shape of a universal dark 
matter profile has been recently proposed by Syer & White 
(1996). In their model, the density profile of a cluster is de- 
termined by the properties of its merging satellites. In par- 
ticular, satellites dense enough to survive the tidal stripping 
of the parent cluster remain self-bound while they sink to 
the cluster centre, dragged by dynamical friction, and are 
therefore responsible for the inner shape of the halo density 
profile. This idea, together with the hierarchical clustering 
model, potentially provides a deeper explanation for the em- 
pirical relation between mass and characteristic density of 
haloes. Since we have reconstructed the merging history of 
all our simulated clusters, it is easy to test the idea directly. 

We first need a definition to identify the matter in the 
dense cores of satellites. For this we take an overdensity 
threshold St, which will always be referred to the mean back- 
ground density at z = 0. This is important, since we need 
to compare the density of satellites identified at different 
redshifts with the density of the final cluster, and we must 
put all haloes on equal footing. For the present analysis we 
took St = 10° ■ Each satellite, identified at any redshift, has 
some fraction of its mass above St- This we shall call the 
mass in satellite cores. We shall say that a satellite is denser 
than another if a larger fraction of its total mass is stored 
in its core, as just defined. By analogy, the mass above the 
threshold S t in the final clusters: Mg t , will be referred to as 
mass in cluster cores. 

We show in Fig. ^| some correlations between different 
quantities derived both from the final clusters and from their 
population of satellites. We used satellites identified at any 
time, and with f c > 0.2. Looking at the figure, we can make 
the following statements. 

(i) Less massive clusters are more centrally concentrated 
(Fig. ^a): the fraction Ms t /M v of cluster mass above the 
fixed overdensity St = 10 5 at z = 0, that is the fraction of 
mass confined in the cluster core, is inversely proportional 
to the total cluster mass M v . 

(ii) Less massive clusters usually form at higher redshift 
(Fig. |§|d): however, the scatter in the relation is fairly big. 
The correlation does not improve if one tries other defini- 
tions of formation time, e.g. changing the fraction of mass 
required in the most massive progenitor from the usual value 
of 50 per cent. 

(iii) The fraction Ms t (sat)/M v of cluster mass provided 
by satellite cores is proportional to the final fraction Mg t /M v 
of mass in the cluster core (Fig. ^s). Therefore, the same 
fraction Mg t (sai)/M v is inversely proportional to the final 
cluster mass (Fig. ^d). 

(iv) The last correlation implies the correlations between 
Ms t (sat)/M v and the scale radius of the profile, r s /R v 
(Fig. ^e). The original correlation between scale radius and 
final cluster mass (as in Tormen et al. 1997, their Fig. 16) is 
reproduced in Fig. ^f. 

Notice that, while the correlation involving the scale radius 
implies a fit to the profile, the others, and in particular that 
of Fig. ^a, is model independent as is directly measured from 
the mass profiles. 
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Figure 8. Correlations between the cluster virial masses M v , the fraction of cluster mass in cluster cores Mi. e s/Mv, the fraction of 
cluster mass which was in satellite cores Mi. e 5(sat)/M v , the cluster formation redshift zj and the scale parameter r s of the cluster density 
profile. M* = (4/3)nR3p b is the characteristic non-linear mass of the model, with p b the mean background density of the universe and 
R t the scale corresponding to a linear matter overdensity 8(R„) = 1.69. Each number labels a cluster. The two straight lines are equal 
weighted-least square fits of the direct (solid) and inverse (dotted) relation. The difference in their slope gives an idea of the goodness of 
the correlation. Definitions of the plotted quantities and details on the figure are given in the text. 



Let us now focus on the correlation of Fig. |8p, and try to 
identify which factors determine it. Is it the mass of satel- 
lites? Their density? Their merging redshift? Their num- 
ber? Or what combination of these? To answer this, we need 
to look in more detail at the characteristics of the satellite 
population. These are shown in Fig. []. We first note, from 
Fig. |^a, that smaller satellites are usually denser, that is a 
higher fraction of their total mass comes from their cores, 
as defined above. This trend is also seen at any fixed time. 
Secondly, Fig. ^|b shows that, on average, smaller clusters 
form from smaller (in physical units) satellites. The error 
bars are ±1<t of the mean values. From these two results we 
may expect that smaller clusters form from denser satellites. 
This is indeed shown in Fig. ^|o. Now, satellites in smaller 
cluster are smaller in absolute terms, but they are actually 
more massive if rescaled to the virial mass of the clusters. 
Therefore, since they are also denser, on average each of 
them brings to the final cluster mass M v a larger contribu- 
tion of mass from satellite cores. This is clearly shown in 
Fig. ^d. The next panel, Fig. ^je, shows that smaller clus- 
ters accrete, on average, a smaller number of satellites. We 
recall that this is not an effect of numerical resolution, as 
each cluster in our sample is resolved by roughly the same 
number of particles, and we consider in all cases satellites 
with at least 8 particles. The average contribution from each 



satellite, given in Fig. ^d, times the total number of satel- 
lites, given in Fig. ^Ja, give exactly the total contribution in 
mass from satellite cores, given in Fig. This shows that 
the higher density of satellites merging onto smaller clusters 
wins over their slightly lower number. 

We finally ask: how relevant to this result is the forma- 
tion redshift of the cluster? If an object forms earlier, it will 
accrete on average smaller, and denser, satellites. One way 
to answer the question is to look at Fig. ^|b: on average, the 
largest progenitor of smaller clusters is assembled earlier. 
However, the distribution of formation redshift predicted by 
the extended P&S formalism is fairly wide, and we noted 
that the correlation in Fig. [sja is not very tight. So we see 
clusters of very different mass which form almost at the same 
time, like clusters labelled 81 and 15, or 40 and 51. How- 
ever, the same pairs of clusters have very different values for 
Ms t /M v , and fairly different values of Ms t (sat)/M v , as well 
as different scale radii for their density profiles (Fig. ^f). A 
different way to address the same question is to measure the 
average redshift at which the cluster final mass was accreted. 
This is just the average value of the mass-weighted identifi- 
cation redshift of the satellite population: < Zid,m v >, with 
m v the satellite virial mass. We can also measure the corre- 
sponding average redshift weighted by the mass in satellite 
cores: < Zidms t >. Both values are shown in Fig. 0f, ver- 
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Figure 9. Correlations between satellite properties and mass of the final cluster. The straight lines are equal weighted-least square fits 
of the direct (solid) and inverse (dashed) relation. The difference in their slope gives an idea of the goodness of the correlation. Each 
number label refers to a cluster. Error bars are ±lcr of the mean. In panel (f), squares are mass weighted by m v , the virial mass of 
satellites. Circles are mass weighted by mg t , the satellite mass above an overdensity threshold St = 10 5 . In this panel, the solid and 
dotted lines refer to the squares, while long dashed and short dashed lines fit the circles. 



sus the final cluster mass M v . The figure shows roughly the 
same trend of Fig. |^b, as expected, but here the scatter in 
the relation is even larger than there. Therefore the redshift 
of formation, or the average redshift of the mergers, is less 
important than the intrinsic density of satellites in produc- 
ing the result of Fig. |^d. 

To summarize, we have seen that, compared to more 
massive clusters, in less massive ones a higher fraction of 
the final mass comes from satellite cores, i.e. is provided by 
denser satellites. These satellites are denser because they are 
smaller, and not necessarily because they are accreted ear- 
lier. Moreover, their cumulative result is due to their density, 
not to their number. We repeated the analysis above using 
different values for the density treshold St, and found that, 
in order to reproduce the result, St must be larger than the 
typical density at the scale radius, but small enough to be 
accessible by the resolution of the simulation. 

To further show that the mass from satellite cores is di- 
rectly responsible for enhancing the central cluster density, 
we compare in Fig. |l^ the density profile of all particles at 
z = (solid curve) to the density profiles, always at z = 0, 
calculated using only the particles from satellite cores. All 
profiles are averaged over the cluster sample. The dotted 
curve refers to satellites merged before the formation red- 
shift Zf, the short-dashed to satellites merged after Zf, and 
the long-dashed to all satellites, regardless to their merging 



time. The figure shows that matter coming from satellite 
cores is more centrally distributed than the average matter 
in the clusters. Among the former, particles merged at earlier 
times have been heated by later mergings, and are less cen- 
trally concentrated. Instead, particles merged more recently 
have a higher central density, but also a larger spreading at 
large radii, possibly because they are less in equilibrium with 
the cluster potential. Quantitatively, the ratio between the 
half-mass radius of the mass from satellite cores and that 
of the whole cluster is rt lm (sat)/rt lm (clu) = 0.45, 1.12, 0.72 
for satellites merged at z > Zf, z < zj and at any z respec- 
tively. No difference is found in the velocity structure of the 
different particle sets. Therefore, particles coming from the 
cores of satellites are more bound to the final system than 
other particles. This confirms the interpretation of the sec- 
ondary infall picture as a conservation of the binding energy 
ranking (Quinn & Zurek 1988; Zaroubi et al. 1996). 

Finally, the ratio Mg t (sat)/M v is only part of the story, 
and other parameters are likely to be important in deter- 
mining the final structure of the cluster profile. For example, 
clusters labelled 23 and 51, which accrete a similar fraction 
of mass in satellite cores, have very different characteris- 
tic density and scale radius (Fig. and f). However, the 
present result makes an important link between the cluster 
dynamical history and its final structure, as it shows that 
the sinking of dense satellites in the centre of their parent 
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Figure 10. Dark matter density profiles at z = 0. The solid curve 
is the profile obtained by using all the particles in the cluster halo. 
The other curves were obtained by using only the particles that 
were in satellite cores. The dotted line is for satellites merged at 
z > Zf, the short dashed line for z < Zf, the long-dashed line 
for all satellites. All profiles are an average over the whole cluster 
sample. Particles from satellite cores are clearly more bound than 
the average particles. Further discussion is given in the text. 

system is a relevant mechanism to giving the cluster its final 
density profile. As far as this mechanism is concerned, less 
massive clusters are more centrally concentrated because a 
larger part of their mass comes from denser satellites. 



cularity e ~ 0.5, avoiding radial orbits as much as almost 
circular ones, and have a mean pericentre at r v ~ 0.38R V . 
More massive haloes tend to fall along more eccentric orbits, 
carrying a smaller specific angular momentum and reaching 
farther into the cluster. 

(ii) The infall of satellites has a very anisotropic distri- 
bution. Its pattern is to a large extent responsible for the 
final shape and orientation of the cluster and of its veloc- 
ity ellipsoid. At redshift z = the major axis of the cluster 
mass distribution is aligned to the major axis of the cluster 
velocity ellipsoid to w 30° on average. The cluster is also 
aligned with the ellipsoid defined by the satellite infall to 
» 30° on average. A similar alignment exists between the 
satellite infall and the final velocity ellipsoid of the cluster. 

Since this infall reflects the whole cluster history after the 
formation redshift at Zf ~ 0.5, the pattern is stable enough 
to provide an explanation for the observed alignment be- 
tween galaxy clusters and the surrounding large-scale struc- 
ture, as recently suggested by the observations of West et 
al. (1995). This result confirms and quantifies the picture, 
often seen in Af-body simulations, of galaxy clusters forming 
from haloes which flow preferentially along filaments. Clus- 
ters shapes and orientations would then reflect the large- 
scale initial density field. The same picture provides a pos- 
sible natural explanation of why massive satellite are ac- 
creted along more eccentric orbits than less massive ones, 
as we found in our orbital analysis. This is because they are 
more likely to reside along dense filaments than less massive 
satellites. Therefore, the filamentary structure acts as a sort 
of focussing rail which drives massive satellites towards the 
cluster. 

(iii) The properties of the population of satellites infalling 
onto the cluster influence also the final shape of its dark mat- 
ter density profile. Lower mass clusters have a larger fraction 
of their total mass supplied by smaller and denser satellites. 
This mechanism is to some extent responsible for the corre- 
lation between the cluster final mass and the characteristic 
overdensity of its profile, observed by Navarro et al. (1996). 



5 SUMMARY AND CONCLUSIONS 

We have analysed the merging of satellites in A'-body sim- 
ulations of a sample of rich galaxy clusters. After identi- 
fication of all cluster progenitors, we selected those haloes 
infalling directly onto the main cluster progenitor. We inves- 
tigated the mass distribution of these satellites, and derived 
the main parameters for their orbits. Finally, we applied 
these results to two problems relevant to cluster formation. 
The first is the relation between the cluster final shape, its 
velocity ellipsoid, and the spatial pattern of the infalling 
satellites, a link that can explain the observed alignment of 
galaxy clusters with the surrounding large-scale structure. 
The second is the correlation between mass and character- 
istic overdensity of dark matter haloes. 

The main conclusions of this paper may be summarized 
as follows. 

(i) The mass function of satellites, i.e. of dark matter 
haloes merging directly onto the main cluster progenitor, 
is in very good agreement with the predictions based on 
Monte Carlo merging trees (Lacey & Cole 1993). Satellites 
fall onto the main progenitor along orbits with average cir- 
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